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Recent groundbreaking experiments studying the effects of spin polarization on pairing in unitary 
Fermi gases encountered mutual qualitative and quantitative discrepancies which seem to be a 
function of the confining geometry. Using novel numerical algorithms we study the solution space for 
a 3-dimensional fully self- consistent formulation of realistic systems with up to 10^ atoms. A study of 
the three types of solutions obtained demonstrates a tendency towards metastability as the confining 
geometry is elongated. One of these solutions, which is consistent with Rice experiments at high trap 
aspect ratio, supports a state strikingly similar to the long sought Fulde-Ferrel-Larkin-Ovchinnikov 
state. Our study helps to resolve the long-standing controversy concerning the discrepancies between 
the findings from two different experimental groups and highlights the versatility of actual-size 
numerical calculations for investigating inhomogeneous fermionic superfluids. 
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Superfluidity in a system of fermionic particles occurs 
when bosonic degrees of freedom emerge and condense 
via pairing of fermions. Understanding the strength of 
this pairing mechanism is closely tied to the search for 
high temperature superconductors. A central issue that 
has animated this quest is: what happens to the pairs 
when the participating species have mismatched Fermi 
surfaces? Such a scenario occurs, for example, in the 
presence of a polarizing field, or when the pairing species 
have unequal numbers or masses. The issue being that, 
when the mismatch of the Fermi surface is large enough, 
a competition between a normal polarized state and the 
superfiuid state would ensue P, [H , potentially giving rise 
to yet unknown or poorly understood exotic superfiuid 
states [3]. Among the interesting theoretical proposals 
for 5- wave pairing is the Fulde-Ferrel-Larkin-Ovchinnikov 
(FFLO) state, a collective term for an inhomogeneous su- 
perfiuid referring either to a Fulde-Ferrell state (FF) 
which supports a super current, or a Larkin-Ovchinnikov 
state (LO) 0) with a spatially modulated order parame- 
ter. Other proposals include breached pairing and p-wave 
symbiotic superfiuids [3]. Primarily due to lack of suffi- 
cient experimental evidence, the issue remains largely un- 
resolved even though it is central to many forms of matter 
such as superconductors, neutron stars and color super- 
fiuids in the quark-gluon plasma [3]. Ultra-cold samples 
of two component degenerate Fermi gases have re- 

energized the debate because of their exquisite controlla- 
bility. 

In this paper, we focus on apparently contradictory re- 
sults on spin-imbalanced unitary Fermi gases from recent 
experiments between two leading groups (7|-[Tq|. In both 
experiments it was observed that, consistent with ear- 
lier predictions [1, 2], the trapped superfiiud responds to 
polarization by phase separating into an inner core with 
negligible polarization surrounded by a polarized outer 
shell. However, the Rice experiments [8, 9] performed in 
cigar-shaped traps with total particle numbers N ~ 10^ 
observed a significant and unexpected deformation of the 



central superfiuid core, indicating a clear violation of the 
local density approximation (LDA). In addition, these 
results also suggest a much higher superfiuid to normal 
(Chandrasekhar-Clogston) transition than the MIT ex- 
periments [7,, J^j in which no deformations were observed. 
The excellent quantitative agreement with theory [ll|, [l2| 
for the MIT experiments conducted at much lower trap 
aspect ratio and with higher particle numbers N ~ 10^, 
hints that there might be unexpected physics at work 
in the Rice experiment. In addition, the concurrence of 
experiments performed in Paris [13] with the MIT ex- 
periments also suggest a crucial role of the trapping ge- 
ometry. This impasse has inspired speculation about the 
possible role of exotic phases such as the FFLO state 
in the observed discrepancies, and stirred much discus- 
sion and debate over the past few years by the cold atom 
community. 

The apparent contradiction between the Rice and 
MIT experiments refiects theoretical difficulties within 
trapped geometries: since the effective chemical potential 
(/i) varies in space, several phases may co-exist within a 
trapped sample. Consequently, despite excitement and 
considerable effort, the theoretical complexity inherent 
within the problem has ensured that most treatments 
have, with few exceptions (l3-[T8|, invoked the LDA 
plol , [20] which is not general enough to capture states 
such as the FFLO. Although an intriguing LDA treat- 
ment which phenomenologically includes a surface energy 
correction has been able to account for the shape of the 
distortions [21], further studies reveal that this model 
is not consistent with a microscopic calculation of the 
surface tension [22]. On the other hand, recent studies 
employing variational techniques in isotropic geometries 
[23, 24] have shown that the region of stability for the 
FFLO state is much larger than originally predicted (25| . 
Until now, a fully self-consistent treatment in anistropic 
geometries with realistic particle numbers, has been well 
out of reach despite its relevance here and in a wide vari- 
ety of other physical systems. To surmount this problem. 
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we developed scalable numerical techniques which take 
full advantage of today's high-performance computing fa- 
cilities running parallel codes over thousands of CPUs. 

We consider a gas of spin-polarized fermionic atoms 
confined to a harmonic trap defined in cylindrical co- 
ordinates {r,(j),z) by V{r,(j),z) = ^(cj^r^ + cj^z^) with 
axial and radial frequencies denoted by {uOz , oo^ ) . Consis- 
tent with Ref. [8, 9] we work at the unitarity limit where 
the 5- wave scattering length between the two spin species 
{as) diverges and within a cigar-shaped trap with aspect 
ratio defined by a = oj±/oJz- This system N = N^-\-Ni 
atoms is described by a Hamiltonian H = f df{Ho -\-Hi) 
with non-interacting {Ho) and interaction {Hi) energy 
densities given by: 



<T=t,4. 

Hi{r) = -U^{f)i;l{r)i;^{r)i;^{f) 



(1) 



where tlJa{'^) and '^J-(r) represents the fermionic field op- 
erators, m the mass and /icr the chemical potential of 
atomic species with spin a. Henceforth, we work in trap 
units for which: m = uJz = ^ = The bare coupling 
constant U is re-normalized through a relationship with 
as by: 1/U = -l/4^a, + (l/F,) E 1/2^^ 112|,126|], with 
Ck = and Vi represents the system volume. H is 

diagonalized through the Bogoliubov-de Gennes (BdG) 
formulation. In particular, our formulation is identical 
to that in Ref. [12]. The superfluid gap (order parame- 
ter) is defined by: A(r) = U{ip^{f)ipi{f)) and the spin 
densities are given by: pa{r) = (V^J {r)tlJa (^))- We find 
it clarifying to express our results in terms of the Fermi 
energy Ep = (3A^)^/^a^/^, and the Thomas- Fermi ra- 
dius along the z-axis Zp = 



\J2Ef for a single species 
ideal Fermi gas of A^/2 particles in a trap with identi- 
cal parameters. In addition, following a convention that 
A^^ > A^;, we define k^^ = -y/2/i^; and the FFLO wave 

number hy qo = k]^ — k'^. 

We solve the BdG equations [l2j using a piece-wise 
linear finite-element basis which yield sparse matrices 
amenable to efficient parallelization and work in a canon- 
ical formalism which fixes N and the total polarization 
P = {N^ — Ni)/N. It has been recently shown that, 
in the particular circumstances of the Rice experiment 
d, evaporative cooling shortens the major axis {z- 
axis) of what should be an ellipsoidal partially polarized 
region, where the condensate forms [23|. Starting from 
an initial ansatz for the gap (A/) imitating this circum- 
stance, the BdG equations are iteratively solved to self- 
consistency using a modified Broyden's method [28]. Our 
calculations reveal that: (1) For large particle numbers 
{N > 10^), we always find a solution similar in struc- 
ture to the LDA solution which has the lowest free en- 
ergy. However starting from an axially shortened initial 
ansatz for the gap, this solution is not accessed by the 
iterative procedure. (2) The most likely solution which 
is consistent with the Rice experiment is a metastable 



state that supports a partially polarized superfluid phase 
strikingly similar to the FFLO phase. This state becomes 
increasingly robust as trapping geometry becomes more 
elongated. (3) Even within a trapped environment, the 
nodes of the order parameter in the FFLO-like phase are 
radially aligned which, with low enough noise, provides 
a measurable, incontrovertible signal within the density 
profiles. 

Superfluidity, a phenomenon of quantum rigidity, ac- 
quires its name from a scenario in which due to energy 
barriers, a condensate gets indefinitely trapped within a 
current-carrying metastable state. The portent for the 
experiments under discussion is that the observed state 
could be a long-lived metastable state. Thus, we take 
the approach of exploring the solution space using an- 
satze constructed with reference to [27] and the phase di- 
agram on the BCS side of the Feshbach resonance j25l , l29| . 
Specifically, we use the LDA solution for the gap (Alda) 
as a base to construct an initial ansatz A/ which is axially 
partitioned into different regions: 

A / ^ I Alda ; \z\ < Zc 

I Alda cos[q{z - z,)] e-(---^)VA^ ; \z\ > z^. 

A I allows us to explore various distorted states. In its 
most general form, one encounters the unpolarized BCS, 
FFLO and normal phases as one traverses along the ax- 
ial direction from the trap center to the edge. The initial 
size of the FFLO region in the ansatz is determined by 
A. When A is too small to accomodate a single wave- 
length of the gap oscillation, i.e., < A < 27r/g', we start 
without an FFLO phase and Zc represents the axial co- 
ordinate of superfluid to normal (S/N) transition. Con- 
versely, an FFLO phase is initially present in the ansatz 
when \ > 2i\ I q. In this case Zc represents the super- 
fluid to FFLO (S/FFLO) transition. Henceforth we refer 
to these initial conditions as Aj~^ and Aj~^^, respec- 
tively, which reflects our nomenclature for the eventual 
solutions as well, i.e., we name the entire solution ac- 
cording to the character of the partially polarized region: 
We have a partially polarized superfluid solution (P-SF) 
when there is an FFLO-like phase present. When the 
partially polarized region is completely normal, we refer 
to the entire solution as a P-N solution. For clarity we 
single out the LDA-like solution which is obtained when 
A/ = Alda, as the SF solution. In both the P-N and 
P-SF solutions, the central unpolarized BCS superfluid 
core is shortened along the z-axis in comparison to the 
LDA-like SF solution. As we shall see, this shortened 
BCS core is manifested in the LDA-violating distortion 
of the density profile of the the minority spin component. 

A broad feature of our results, which directly informs 
on the question of metastability, is the observation of a 
barrier between the shortened states (either P-N or P- 
SF) and the SF solution. For small atom numbers, this 
barrier is absent, the converged solution is unique, in- 
dependent of the initial ansatz we take, and we see a 
dramatic departure from the LDA prediction due to sig- 
nificant finite-size effect. However with increasing A", the 
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We ascribe the consistent convergence to a shortened 
state as due to the emergence of energy barriers sepa- 
rating the P-SF and P-N states from the SF state with 
increasing or a in tandem with Ep. In Fig. [TJa) we 
ihustrate the dramatic differences in the superfluid gap 
for the various solutions encountered. Apart from the 
emerging energy barriers, another important result with 
regard to metastability is the decrease in the relative en- 
ergetic separation of all the states, P-SF, P-N and SF, as 
a is increased. Taken together, these observations sug- 
gest that the relaxation of the physical system from any 
of the shortened states to the SF state, which is the low- 
est in energy, becomes less favorable as a is increased, a 
deduction which is borne out by the discrepancies of the 
Rice and MIT experiments. 

For a given value of Zc, the energy of the P-SF solution 
is consistently lower than the P-N solution. Furthermore, 
recent results suggest that the inclusion of fluctuations, 
neglected in mean- field formulations, should make the P- 
SF state even more stable (23|. Thus, we expect that if 
the system converges to a shortened state, it will choose 
the P-SF state. A natural question to ask is: how will 
the FFLO phase manifests itself? 



Figure 1: (Color online) (a) Axial profiles of the gap (in units 
of ^f) showing the P-N (blue dashed line), P-SF (red solid 
line) and SF (green dotted line) states. The LDA solution 
(not shown) almost completely overlaps with the SF result. 
The free energies per particle are: 0.67(0)£^f, 0.65(8)£^f, 
0.65(5)^F and 0.64(4)^f for the P-N, P-SF, SF and LDA 
states, respectively, (b) Local polarization p(r) within the 
partially polarized region of the P-SF (red solid line) and P- 
N(blue dashed line) solutions, (c) An r-z plot of the nor- 
malized density difference bp — [p^ — pi)/pF of the partially 
polarized region of the P-SF state (pp = yj {2EfY /Qi^'^). All 
the results shown in this paper are obtained at a small tem- 
perature T = ^mEp/kB, and with N = 50000, a = 50, and 
P = 0.3. 
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Figure 2: (color online) Plots showing the doubly integrated 
axial spin density 8pid{z) for, from top down, the SF, P-N 
and P-SF states shown in Fig. [TJa). 



axial S/N or S/FFLO transition point is pinned near its 
initial value Zc and we obtain different solutions by start- 
ing from different initial ansatze. Starting from Aj or 
Aj ~^ we always converge to a shortened state in a man- 
ner which is only sensitive to our choice oi q. In other 
words, we do see a transition between the P-SF and P-N 
states which is very sensitive to q and largely insensitive 
to A; both of which are set in the initial condition A/. It 
works as follows. When q is less than a critical value qc^ 
the oscillations in the ansatz A/ are amplified and the 
solution flows to a P-SF state no matter the size of A. 
Conversely when q > qc, the oscillations are damped and 
A/ always converges to a P-N state. A similar resonance 
behavior has also recently been observed in studies of 
the S/N boundary while tuning across the BEC-BCS 
crossover [22], in which case calculations were performed 
without the radial confinement. It is possible that this 
phenomenon might be exploited to engineer the realiza- 
tion of the P-SF state. 



In Fig. [TJb) we contrast the appearance of local polar- 
ization p{r) = {p^—pi) I {p^^p\^) in the partially polarized 
regions of the P-SF and P-N states. We note that in [23|, 
the strong oscillations displayed in p{r) were observed to 
survive the effects of fluctuations. One pleasant surprise 
of our results was the radial alignment of the nodes of 
the FFLO phase shown in Fig. Hl^c). A fact which is not 
a priori obvious, and very promising for the prospects 
of detection within the 3D system under discussion here, 
because it implies that the FFLO phase could yield a 
measurable signal in the density profile. Auspiciously, it 
also suggests that when an array of ID tubes, such as 
are being used in current experiments [30], are coupled 
to yield a quasi-3D confinement, the FFLO nodes at each 
tube are likely to align to yield a measurable signal. To 
make sure that the radial alignment of the nodes is not 
a numerical artefact, we have used initial ansatze where 
the nodes are intentionally misaligned along the radial 
direction. Our code always converges to states with the 
nodes aligned. A comparison of the plots in Fig. [2] con- 
firms that the presence of an FFLO phase would indeed 
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provide a smoking gun signal in doubly integrated axial 
spin density dpid = f f dxdy{p^ — p^). In the close-up 
we observe that the signal of the FFLO region is not as 
strong as that in Fig. [TJc) because of contributions from 
the fully polarized shell encasing it. Quantitatively, it 
indicates that a lower bound of the signal to noise ratio 
~ 6.5 is required to observe at least half of the FFLO 
phase. 

A casual comparison of all column density profiles in 
Fig. [3] rules out the observation in Rice experiment of the 
SF state, which is consistent with the LDA and, within 
the BdG formulation, has the lowest free energy. How- 
ever, due to the noise on the experimental data, it is not 
clear which of the shortened states (P-SF or P-N) has 
been observed. To produce the noise with similar char- 
acteristics as the experiment, we added white noise with 
standard deviation which is a similar fraction of the aver- 
age value of the column density p^ dx in the plotted 
window. Theoretically, since it has the lower energy and 
since the transition between the FFLO phase and the nor- 
mal phase is continuous, one expects that, between the 
two shortened states, the P-SF solution will be favored. 

In conclusion, we have repeatedly solved the BdG 
equations in a cigar-shaped trap using initial conditions 
which imitate the condensate nucleation process (27| . 
The iterative solution chooses between two stationary 
points, which are not necessarily the global free energy 
minimum, but each of which features density profiles 
strikingly similar to experimental observations at Rice. 
The solution which possesses the lower energy of these 
two contains an FFLO-like phase which leaves an acces- 
sible signal in 5pid. Coupled with recent results which 
suggest the unexpected stability of the FFLO in 3D [24], 
our observations raise the interesting question of whether 



the FFLO state has already been realized in the Rice ex- 
periment. Since the Hartree interactions are excluded 
from the BdG formulation of unitary gases [12], we do 
not address the position of the Clogston limit. Never- 
theless we note in passing that the P-SF solution has the 
capacity to absorb polarizations and, if undetected, could 
conceal the existence of a partially polarized region. Fi- 
nally we remark that our work is important for another 
reason, as far as we know ours are by far the largest 
calculations of their kind and herald the arrival of an 
important tool for investigations of finite fermionic sys- 
tems such as occur in atomic traps or in nuclear physics; 
where predictions of ideal models such as the FFLO pro- 
posal could be significantly modified by confinement and 
finite-size effects. 

Note added After our work was completed, the Rice ex- 
perimental group verified the suggestion made in Ref. |23| 
that the LDA-violating deformations observed in their 
experiment are a result of depolarization of the super- 
fluid core by evaporation occurring mainly at the axial 
center of the trap fsij. They found that these deformed 
states are very stable, in agreement with our calculations. 
The metastability of these states suggests the possiblity 
to directly engineer an FFLO state in an elongated trap. 
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Figure 3: (color online) Each display shows the column densities (rescaled to have aspect ratio 5 for clarity) / p-^dx, J p^dx, 
J{p^ — pi)dx and the axial spin density 6pid, respectively. The states represented are (a) the SF, (b) P-SF and (c) P-N states 
illustrated in Fig.IHa). In (d) we plot the Rice experimental results for N ^ 260000, P ^ 0.35, a = 45.23 and T < OmEp/kB- 
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